module cubed_sphere_domain_mod

  use mpi
  use flogger
  use kinds_mod
  use mesh_const_mod
  use mesh_math_mod
  use cubed_sphere_math_mod

  implicit none

  private

  public cubed_sphere_domain_type

  integer, public, parameter :: equiangular = 1

  type ghost_point_type
    integer :: ngp = 0 ! Number of ghost points on this cell
    
    integer    , allocatable, dimension(:) :: idom ! Source domain index
    integer    , allocatable, dimension(:) :: ps   ! Source point index
    integer    , allocatable, dimension(:) :: is   ! Source cell index
    integer    , allocatable, dimension(:) :: js   ! Source cell index
    real   (r8), allocatable, dimension(:) :: wgt  ! Weight to take into account of point resides on domain edges
    
    real   (r8), allocatable, dimension(:) :: x
    real   (r8), allocatable, dimension(:) :: y
    real   (r8), allocatable, dimension(:) :: lon
    real   (r8), allocatable, dimension(:) :: lat
           
    real   (r8), allocatable, dimension(:,:) :: r
    real   (r8), allocatable, dimension(:,:) :: z
    real   (r8), allocatable, dimension(:,:) :: eta
    real   (r8), allocatable, dimension(:,:) :: xi
    
    ! Dimension columns:
    ! 1. row
    ! 2. colume
    ! 3. point
    ! 4. z-axis
    !                                   1 2 3 4
    real   (r8), allocatable, dimension(:,:,:  ) :: Ah    ! Cube-to-sphere transformation matrix (Gh = Ah^T Ah)
    real   (r8), allocatable, dimension(:,:,:  ) :: iAh   ! Sphere-to-cube transformation matrix
    real   (r8), allocatable, dimension(:,:,:  ) :: Gh    ! Metric tensor (contravariant-to-covariant vector transformation matrix)
    real   (r8), allocatable, dimension(:,:,:  ) :: iGh   ! Inversed metric tensor (covariant-to-contravariant vector transformation matrix)
    real   (r8), allocatable, dimension(    :  ) :: Jh    ! Jacobian (det(Gv))
    
    real   (r8), allocatable, dimension(    :,:) :: Jv    ! Jacobian (det(Gv))
    
    real   (r8), allocatable, dimension(    :,:) :: Jab   ! Jacobian (det(G)), Jh * Jv
    contains
    procedure :: clear => ghost_clear
    final     :: ghost_final
  end type ghost_point_type

  type cubed_sphere_ngb_domain_type
    type(cubed_sphere_domain_type), pointer :: ptr => null()
    integer :: orient = 0
    integer :: cell_type
    integer :: code = 0
    integer :: ex_depth = 0 ! Depth extension
  end type cubed_sphere_ngb_domain_type

  type cubed_sphere_domain_type
    logical  :: initialized = .false.
    logical  :: active      = .false.
    integer  :: proj_type
    integer  :: domain_type = global_domain
    integer  :: pid  = MPI_PROC_NULL
    integer  :: id   = 0
    integer  :: hdim = 0                  ! mesh dimension, choose from 2 or 3
    integer  :: nc   = 0                  ! Resolution parameter
    integer  :: nx   = 0                  ! Total number of cells along x-axis
    integer  :: ny   = 0                  ! Total number of cells along y-axis
    integer  :: nz   = 1                  ! Number of vertical levels
    integer  :: ndx  = 0                  ! Number of domain cells along x-axis
    integer  :: ndy  = 0                  ! Number of domain cells along y-axis
    integer  :: nbx  = 0                  ! Number of block cells along x-axis
    integer  :: nby  = 0                  ! Number of block cells along y-axis
    integer  :: hw0  = 0                  ! Base halo width
    integer  :: hw(4)= 0                  ! Halo width along each side, 1 left, 2 right, 3 bottom, 4 top
    integer  :: ex(8)= 0                  ! Extension along each direction
    integer  :: rw   = 0                  ! Horizontal Reconstruction width
    integer  :: rwv  = 0                  ! Vertical Reconstruction width
    real(r8) :: x0   = 0._r8              ! Cell center start x
    real(r8) :: y0   = 0._r8              ! Cell center start y
    real(r8) :: xds  = 0._r8              ! Start x on process
    real(r8) :: xde  = 0._r8              ! End x on process
    real(r8) :: yds  = 0._r8              ! Start y on process
    real(r8) :: yde  = 0._r8              ! End y on process
    real(r8) :: xmin = 0._r8              ! Start x on panel
    real(r8) :: xmax = 0._r8              ! End x on panel
    real(r8) :: ymin = 0._r8              ! Start y on panel
    real(r8) :: ymax = 0._r8              ! End y on panel
    real(r8) :: zmin = 0._r8              ! min z
    real(r8) :: zmax = 0._r8              ! max z
    real(r8) :: etamin = 0._r8            ! min eta
    real(r8) :: etamax = 0._r8            ! max eta
    real(r8) :: dx   = 0._r8              ! Computational grid interval along x-axis
    real(r8) :: dy   = 0._r8              ! Computational grid interval along y-axis
    real(r8) :: deta = 0._r8              ! Computational grid interval along z-axis
    ! Cell index ranges
    integer :: ids  = 0, ide  = 0         ! Cell start and end indices along x-axis
    integer :: ims  = 0, ime  = 0         ! Cell with halo start and end indices along x-axis
    integer :: jds  = 0, jde  = 0         ! Cell start and end indices along y-axis
    integer :: jms  = 0, jme  = 0         ! Cell with halo start and end indices along y-axis
    integer :: kms  = 0, kme  = 0         ! Cell with halo start and end indices along y-axis
    integer :: kds  = 0, kde  = 0         ! Cell with halo start and end indices along y-axis
    integer :: ipds = 0, ipde = 0         ! Index for in-domain cells along x-axis
    integer :: jpds = 0, jpde = 0         ! Index for in-domain cells along y-axis
    integer :: irs  = 0, ire  = 0         ! Index for the cells required by reconstruction along x-axis
    integer :: jrs  = 0, jre  = 0         ! Index for the cells required by reconstruction along y-axis
    integer :: igs  = 0, ige  = 0         ! Index for ghost cells along x-axis
    integer :: jgs  = 0, jge  = 0         ! Index for ghost cells along y-axis
    integer :: its  = 0, ite  = 0         ! Index for reconstruction target cell along x-axis
    integer :: jts  = 0, jte  = 0         ! Index for reconstruction target cell along y-axis
    ! Point type and their indicators
    integer :: npt                        ! Number of points in each cell (e.g., cell center, vertex)
    integer :: nv                         ! Number of vertices in cell
    integer :: ncq                        ! Number of quadrature points in cell
    integer :: nvq                        ! Number of quadrature points in cell volume
    integer :: neq                        ! Number of quadrature points along each horizontal edge
    integer :: nfq                        ! Number of quadrature points on each cell face
    integer :: nfqv                       ! Number of quadrature points along each vertical edge
    integer :: pc   = 1                   ! Point index for cell center
    integer :: pvs                        ! Point start index for vertices
    integer :: pve                        ! Point end index for vertices
    integer :: pes(6)                     ! Point start index for quadrature points on edges
    integer :: pee(6)                     ! Point end index for quadrature points on edges
    integer :: pqs                        ! Point start index for quadrature points in cell
    integer :: pqe                        ! Point end index for quadrature points in cell
    integer :: pgs                        ! Point start index for ghost points in cell
    integer :: pge                        ! Point maximum end index for ghost points in cell
    ! Quadrature points
    real(16), allocatable :: xeq(:)       ! Quadrature locations
    real(r8), allocatable :: weq1d(:)     ! Quadrature weights for 1D
    real(r8), allocatable :: weq2d(:)     ! Quadrature weights for 2D
    real(r8), allocatable :: weq3d(:)     ! Quadrature weights for 3D
    real(r8), allocatable :: weq_v(:)     ! Vertical quadrature weights
    real(r8), allocatable :: weq_edge(:)  ! Quadrature weights for edge
    real(r8), allocatable :: weq_cell(:)  ! Quadrature weights for cell
    ! Dimension columns:
    ! 1. row
    ! 2. colume
    ! 3. colume 2, for Christoffel Symbol only
    ! 4. point
    ! 5. x-axis
    ! 6. y-axis
    ! 7. z-axis
    !                                1 2 3 4 5 6 7
    real(r8), allocatable, dimension(      :,:,:,:) :: r     ! radius vector length
    real(r8), allocatable, dimension(      :,:,:  ) :: x     ! Point x coordinate
    real(r8), allocatable, dimension(      :,:,:  ) :: y     ! Point y coordinate
    real(r8), allocatable, dimension(      :,:,:,:) :: z     ! Point z coordinate
    real(r8), allocatable, dimension(      :,    :) :: eta   ! Point eta coordinate, vertical coordinate in computational space
    real(r8), allocatable, dimension(      :,    :) :: xi    ! Point xi coordinate, vertical height in computational space
    real(r8), allocatable, dimension(      :,:,:  ) :: lon   ! Point longitude
    real(r8), allocatable, dimension(      :,:,:  ) :: lat   ! Point latitude
    
    real(r8), allocatable, dimension(:,:,  :,:,:  ) :: Ah    ! Cube-to-sphere transformation matrix (Gh = Ah^T Ah)
    real(r8), allocatable, dimension(:,:,  :,:,:  ) :: iAh   ! Sphere-to-cube transformation matrix
    real(r8), allocatable, dimension(:,:,  :,:,:  ) :: Gh    ! Metric tensor (contravariant-to-covariant vector transformation matrix)
    real(r8), allocatable, dimension(:,:,  :,:,:  ) :: iGh   ! Inversed metric tensor (covariant-to-contravariant vector transformation matrix)
    real(r8), allocatable, dimension(      :,:,:  ) :: Jh    ! Jacobian (det(Gh))
    real(r8), allocatable, dimension(:,:,:,:,:,:  ) :: CS    ! Christoffel Symbol 2nd type
    
    real(r8), allocatable, dimension(:,:,  :,:,:,:) :: iGv   ! Inversed Gv
    real(r8), allocatable, dimension(      :,:,:,:) :: Jv    ! Jacobian (det(Gv))
    
    real(r8), allocatable, dimension(      :,:,:,:) :: Jab   ! Jacobian (det(G)), Jh * Jv
    real(r8), allocatable, dimension(:,:,  :,:,:,:) :: iG    ! Inversed G
    
    real(r8), allocatable, dimension(:,    :,:,:,:) :: omg   ! planet self-roation angle velocity

    ! Location indicators
    logical :: touch_edges(4)  = .false.
    integer, allocatable, dimension(:,:) :: cell_type

    ! Ghost points
    ! 1. x-axis
    ! 2. y-axis
    !                                              1 2
    type(ghost_point_type), allocatable, dimension(:,:) :: ghost
    integer :: n_gst_cells
    integer, dimension(:,:), allocatable :: gst_ij

    integer :: n_halo_cells
    integer, dimension(:,:), allocatable :: halo_ij

    integer :: max_ngp = 0
    ! Neighbor domains
    integer :: ngb_code = 0
    type(cubed_sphere_ngb_domain_type) ngb(4)
    ! Coordinate transformation procedures
    procedure(locate_domain_interface     ), pointer, nopass :: locate_domain
    procedure(cube_to_sphere_interface    ), pointer, nopass :: cube_to_sphere
    procedure(sphere_to_cube_interface    ), pointer, nopass :: sphere_to_cube
    procedure(horizontal_metrics_interface), pointer, nopass :: horizontal_metrics
    procedure(vertical_metrics_interface  ), pointer, nopass :: vertical_metrics
  contains
    procedure          :: connect               => cubed_sphere_domain_connect
    procedure          :: init                  => cubed_sphere_domain_init
    procedure, private :: init_equiangular      => cubed_sphere_domain_init_equiangular
    procedure          :: set_hor_metrics       => cubed_sphere_domain_set_horizontal_metrics
    procedure          :: set_ver_metrics       => cubed_sphere_domain_set_vertical_metrics
    procedure          :: set_ghost_ver_metrics => cubed_sphere_domain_set_ghost_vertical_metrics
    procedure, private :: count_ghost_point     => cubed_sphere_domain_count_ghost_point
    procedure          :: count_ghost_points    => cubed_sphere_domain_count_ghost_points
    procedure          :: allocate_ghost_points => cubed_sphere_domain_allocate_ghost_points
    procedure, private :: add_ghost_point       => cubed_sphere_domain_add_ghost_point
    procedure          :: add_ghost_points      => cubed_sphere_domain_add_ghost_points
    procedure          :: set_ghost_ij          => cubed_sphere_domain_set_ghost_ij
    procedure          :: edge_quad             => cubed_sphere_domain_edge_quad
    procedure          :: edge_quad_ver         => cubed_sphere_domain_edge_quad_ver
    procedure          :: cell_quad             => cubed_sphere_domain_cell_quad
    procedure          :: is_xy_in_domain       => cubed_sphere_domain_is_xy_in_domain
    procedure          :: clear                 => cubed_sphere_domain_clear
    final :: cubed_sphere_domain_final
  end type cubed_sphere_domain_type

contains

  subroutine cubed_sphere_domain_connect(this, id, left_domain, right_domain, top_domain, bottom_domain)

    class(cubed_sphere_domain_type), intent(inout) :: this
    integer, intent(in) :: id
    type(cubed_sphere_domain_type), intent(in), target :: left_domain
    type(cubed_sphere_domain_type), intent(in), target :: right_domain
    type(cubed_sphere_domain_type), intent(in), target :: top_domain
    type(cubed_sphere_domain_type), intent(in), target :: bottom_domain

    this%id = id
    this%ngb(1)%ptr =>   left_domain
    this%ngb(2)%ptr =>  right_domain
    this%ngb(3)%ptr => bottom_domain
    this%ngb(4)%ptr =>    top_domain

    ! Domain topology
    !        _____
    !       |     |
    !       |  5  |
    !  _____|_____|___________
    ! |     |     |     |     |
    ! |  4  |  1  |  2  |  3  |
    ! |_____|_____|_____|_____|
    !       |     |
    !       |  6  |
    !       |_____|
    !
    select case (id)
    case (1)
      this%ngb(1)%orient = right ; this%ngb(1)%code = 14
      this%ngb(2)%orient = left  ; this%ngb(2)%code = 12
      this%ngb(3)%orient = top   ; this%ngb(3)%code = 16
      this%ngb(4)%orient = bottom; this%ngb(4)%code = 15
    case (2)
      this%ngb(1)%orient = right ; this%ngb(1)%code = 21
      this%ngb(2)%orient = left  ; this%ngb(2)%code = 23
      this%ngb(3)%orient = right ; this%ngb(3)%code = 26
      this%ngb(4)%orient = right ; this%ngb(4)%code = 25
    case (3)
      this%ngb(1)%orient = right ; this%ngb(1)%code = 32
      this%ngb(2)%orient = left  ; this%ngb(2)%code = 34
      this%ngb(3)%orient = bottom; this%ngb(3)%code = 36
      this%ngb(4)%orient = top   ; this%ngb(4)%code = 35
    case (4)
      this%ngb(1)%orient = right ; this%ngb(1)%code = 43
      this%ngb(2)%orient = left  ; this%ngb(2)%code = 41
      this%ngb(3)%orient = left  ; this%ngb(3)%code = 46
      this%ngb(4)%orient = left  ; this%ngb(4)%code = 45
    case (5)
      this%ngb(1)%orient = top   ; this%ngb(1)%code = 54
      this%ngb(2)%orient = top   ; this%ngb(2)%code = 52
      this%ngb(3)%orient = top   ; this%ngb(3)%code = 51
      this%ngb(4)%orient = top   ; this%ngb(4)%code = 53
    case (6)
      this%ngb(1)%orient = bottom; this%ngb(1)%code = 64
      this%ngb(2)%orient = bottom; this%ngb(2)%code = 62
      this%ngb(3)%orient = bottom; this%ngb(3)%code = 63
      this%ngb(4)%orient = bottom; this%ngb(4)%code = 61
    end select
    this%ngb(1)%cell_type =   left_halo_cell
    this%ngb(2)%cell_type =  right_halo_cell
    this%ngb(3)%cell_type = bottom_halo_cell
    this%ngb(4)%cell_type =    top_halo_cell

  end subroutine cubed_sphere_domain_connect

  recursive subroutine cubed_sphere_domain_init(this, proj_type, hdim, vdim, nc, ndx, ndy, hw, hws, rw, rwv, &
                                                neq, xeq, weq, ids, jds, nbx, nby, nz, active,  &
                                                ngb_code, pid, cart_dims, cart_coords)

    class(cubed_sphere_domain_type), intent(inout) :: this
    integer , intent(in)           :: proj_type      ! Projection type (only equiangular for now)
    integer , intent(in), optional :: hdim           ! horizontal reconstruction dimension, choose from 2 or 3
    integer , intent(in), optional :: vdim           ! vertical reconstruction dimension, choose from 1 to 3
    integer , intent(in), optional :: nc             ! Resolution parameter
    integer , intent(in), optional :: ndx            ! Number of domain cells along x-axis
    integer , intent(in), optional :: ndy            ! Number of domain cells along y-axis
    integer , intent(in), optional :: hw             ! Halo width
    integer , intent(in), optional :: hws(4)         ! Halo width array
    integer , intent(in), optional :: rw             ! Horizontal Reconstruction width
    integer , intent(in), optional :: rwv            ! Vertical Reconstruction width
    integer , intent(in), optional :: neq            ! Number of quadrature points along cell edge
    real(16), intent(in), optional :: xeq(:)         ! Positions of quadrature points along cell edge
    real(r8), intent(in), optional :: weq(:)         ! Weights of quadrature points along cell edge
    integer , intent(in), optional :: ids            ! Start index along x axis
    integer , intent(in), optional :: jds            ! Start index along y axis
    integer , intent(in), optional :: nbx            ! Number of grids along x axis in block
    integer , intent(in), optional :: nby            ! Number of grids along y axis in block
    integer , intent(in), optional :: nz             ! Number of vertical levels
    logical , intent(in), optional :: active         ! Flag for indicating this is the main domain in the process
    integer , intent(in), optional :: ngb_code       ! Neighbor code
    integer , intent(in), optional :: pid            ! Process ID
    integer , intent(in), optional :: cart_dims(2)   ! Process Cartesian decomposition dimensions
    integer , intent(in), optional :: cart_coords(2) ! Process Cartesian decomposition coordinates
    
    integer is, ie, js, je, depth, ngb_ids, ngb_jds, ngb_nx, ngb_ny, ngb_hw(4)
    integer ex_depth, ex_width, max_ex_depth

    call this%clear()

    max_ex_depth = 6
    
    if(present(hdim))then
      this%hdim = hdim
    else
      this%hdim = 2
    endif
    
    if (present(nc)) then
      this%nc  = nc
      this%ndx = nc
      this%ndy = nc
    else if (present(ndx) .and. present(ndy)) then
      this%domain_type = regional_domain
      this%nc  = 0
      this%ndx = ndx
      this%ndy = ndy
    end if
    if (present(nz )) this%nz  = nz
    if (present(pid)) this%pid = pid

    if (present(ngb_code)) this%ngb_code = ngb_code
    
    if(this%hdim==2)then
      this%pvs = 2
      this%pve = 5
    elseif(this%hdim==3)then
      this%pvs = 2
      this%pve = 9
    endif
    
    ! Set indices for each type of point.
    if (present(neq)) then
      this%neq = neq
    else
      this%neq = 1
    end if

    if(this%hdim==2)then
      this%neq = neq
      this%nfq = this%neq
      this%ncq = this%neq**2
      this%nvq = this%neq**2
      this%nv  = 4
    elseif(this%hdim==3)then
      this%neq = neq
      this%nfq = this%neq**2
      this%ncq = this%neq**2
      this%nvq = this%neq**3
      this%nv  = 8
    endif
    
    if( vdim == 1 )then
      this%nfqv = 1
    else
      this%nfqv = this%ncq
    endif
    
    this%max_ngp = 4 * this%nvq
    
    this%pes(3) = this%pve + 1 + 0*this%nfq
    this%pee(3) = this%pve + 1 + 1*this%nfq - 1
    this%pes(2) = this%pve + 1 + 1*this%nfq
    this%pee(2) = this%pve + 1 + 2*this%nfq - 1
    this%pes(4) = this%pve + 1 + 2*this%nfq
    this%pee(4) = this%pve + 1 + 3*this%nfq - 1
    this%pes(1) = this%pve + 1 + 3*this%nfq
    this%pee(1) = this%pve + 1 + 4*this%nfq - 1
    
    ! Down Face Points
    this%pes(5) = this%pee(1) + 1
    this%pee(5) = this%pes(5) + this%nfqv - 1
    ! Up Face Points
    this%pes(6) = this%pee(5) + 1
    this%pee(6) = this%pes(6) + this%nfqv - 1
    
    this%pqs = this%pee(6) + 1
    this%pqe = this%pqs + this%nvq - 1
    
    this%npt = this%pqe
    
    ! Ghost point is added later.
    this%pgs = this%pqe + 1
    this%pge = this%pgs + this%max_ngp - 1
    
    allocate(this%xeq  (neq   ))
    allocate(this%weq1d(neq   ))
    allocate(this%weq2d(neq**2))
    allocate(this%weq3d(neq**3))
    allocate(this%weq_v(this%nfqv))

    if( hdim==2 )then
      allocate(this%weq_edge(neq**1))
      allocate(this%weq_cell(neq**2))
    elseif( hdim==3 )then
      allocate(this%weq_edge(neq**2))
      allocate(this%weq_cell(neq**3))
    end if

    if (present(neq) .and. .not. (present(xeq) .and. present(weq))) then
      ! Generate Gaussian-Legendre quadrature points for users.
      call gaussian_legendre(neq, this%xeq, this%weq1d)
      this%xeq   = (this%xeq + 1) * 0.5_16
      this%weq1d = this%weq1d * 0.5_r8
      this%weq2d = kronecker_product(neq,this%weq1d,neq,this%weq1d)
      this%weq3d = kronecker_product(neq**2,this%weq2d,neq,this%weq1d)
    else if (present(neq) .and. present(xeq) .and. present(weq)) then
      this%xeq   = xeq
      this%weq1d = weq
      this%weq2d = kronecker_product(neq,this%weq1d,neq,this%weq1d)
      this%weq3d = kronecker_product(neq**2,this%weq2d,neq,this%weq1d)
    end if

    if( hdim==2 )then
      this%weq_edge = this%weq1d
      this%weq_cell = this%weq2d
    elseif( hdim==3 )then
      this%weq_edge = this%weq2d
      this%weq_cell = this%weq3d
    end if
    
    if( vdim == 1 )then
      this%weq_v = 1
    else
      this%weq_v = this%weq2d
    endif

    this%active      = merge(active, .false., present(active))
    this%proj_type   = proj_type
    this%domain_type = merge(regional_domain, global_domain, this%nc == 0)
    
    if (present(hw)) then
      this%hw0 = hw
      this%hw = hw
    else if (present(hws)) then
      this%hw = hws
    end if
    
    this%nx = this%ndx + this%hw(1) + this%hw(2)
    this%ny = this%ndy + this%hw(3) + this%hw(4)

    if (present(rw)) then
      this%rw = rw
    else
      this%rw = this%hw(1)
    end if
    depth = 2 * this%rw + 1
    
    if(present(rwv))then
      this%rwv = rwv
    else
      this%rwv = 0
    end if

    select case (proj_type)
    case (equiangular)
      this%locate_domain      => equiangular_locate_domain
      this%cube_to_sphere     => equiangular_cube_to_sphere
      this%sphere_to_cube     => equiangular_sphere_to_cube
      this%horizontal_metrics => equiangular_metrics
      this%vertical_metrics   => vertical_metrics
    end select

    if (present(ids) .and. present(nbx)) then
      this%ids = ids; this%ide = ids + nbx - 1; this%nbx = nbx
    else
      this%ids = 1  ; this%ide = nc           ; this%nbx = nc
    end if
    if (present(jds) .and. present(nby)) then
      this%jds = jds; this%jde = jds + nby - 1; this%nby = nby
    else
      this%jds = 1  ; this%jde = nc           ; this%nby = nc
    end if

    this%touch_edges(1) = this%ids == 1
    this%touch_edges(2) = this%ide == this%nc
    this%touch_edges(3) = this%jds == 1
    this%touch_edges(4) = this%jde == this%nc

    ex_depth = this%rw
    ex_width = 4

    if (this%active) then
      if (all(this%touch_edges)) then
        this%ex = [ex_depth,ex_depth,ex_depth,ex_depth,0,0,0,0]
      ! Left-bottom corner
      else if (this%touch_edges(1) .and. .not. this%touch_edges(2) .and. this%touch_edges(3) .and. .not. this%touch_edges(4)) then
        this%hw = [hw,hw+ex_width,hw,hw+ex_width]
        this%ex = [ex_depth,0,ex_depth,0,0,ex_depth,0,ex_depth]
        if (present(cart_dims) .and. present(cart_coords)) then
          if (cart_coords(1) + 1 < cart_dims(1) - 1) this%ex(right_bottom) = 0
          if (cart_coords(2) + 1 < cart_dims(2) - 1) this%ex(left_top) = 0
        end if
      ! Bottom inner
      else if (.not. this%touch_edges(1) .and. .not. this%touch_edges(2) .and. this%touch_edges(3) .and. .not. this%touch_edges(4)) then
        this%hw = [hw+ex_width,hw+ex_width,hw,hw]
        this%ex = [0,0,0,0,0,0,0,0]
      ! Right-bottom corner
      else if (.not. this%touch_edges(1) .and. this%touch_edges(2) .and. this%touch_edges(3) .and. .not. this%touch_edges(4)) then
        this%hw = [hw+ex_width,hw,hw,hw+ex_width]
        this%ex = [0,ex_depth,ex_depth,0,ex_depth,0,ex_depth,0]
        if (present(cart_dims) .and. present(cart_coords)) then
          if (cart_coords(1) - 1 > 0) this%ex(left_bottom) = 0
          if (cart_coords(2) + 1 < cart_dims(2) - 1) this%ex(right_top) = 0
        end if
      ! Right inner
      else if (.not. this%touch_edges(1) .and. this%touch_edges(2) .and. .not. this%touch_edges(3) .and. .not. this%touch_edges(4)) then
        this%hw = [hw,hw,hw+ex_width,hw+ex_width]
        this%ex = [0,0,0,0,0,0,0,0]
      ! Right-top corner
      else if (.not. this%touch_edges(1) .and. this%touch_edges(2) .and. .not. this%touch_edges(3) .and. this%touch_edges(4)) then
        this%hw = [hw+ex_width,hw,hw+ex_width,hw]
        this%ex = [0,ex_depth,0,ex_depth,0,ex_depth,0,ex_depth]
        if (present(cart_dims) .and. present(cart_coords)) then
          if (cart_coords(1) - 1 > 0) this%ex(left_top) = 0
          if (cart_coords(2) - 1 > 0) this%ex(right_bottom) = 0
        end if
      ! Top inner
      else if (.not. this%touch_edges(1) .and. .not. this%touch_edges(2) .and. .not. this%touch_edges(3) .and. this%touch_edges(4)) then
        this%hw = [hw+ex_width,hw+ex_width,hw,hw]
        this%ex = [0,0,0,0,0,0,0,0]
      ! Left-top corner
      else if (this%touch_edges(1) .and. .not. this%touch_edges(2) .and. .not. this%touch_edges(3) .and. this%touch_edges(4)) then
        this%hw = [hw,hw+ex_width,hw+ex_width,hw]
        this%ex = [ex_depth,0,0,ex_depth,ex_depth,0,ex_depth,0]
        if (present(cart_dims) .and. present(cart_coords)) then
          if (cart_coords(1) + 1 < cart_dims(1) - 1) this%ex(right_top) = 0
          if (cart_coords(2) - 1 > 0) this%ex(left_bottom) = 0
        end if
      ! Left inner
      else if (this%touch_edges(1) .and. .not. this%touch_edges(2) .and. .not. this%touch_edges(3) .and. .not. this%touch_edges(4)) then
        this%hw = [hw,hw,hw+ex_width,hw+ex_width]
        this%ex = [0,0,0,0,0,0,0,0]
      end if
    end if

    this%ims = this%ids - this%hw(1); this%ime = this%ide + this%hw(2)
    this%jms = this%jds - this%hw(3); this%jme = this%jde + this%hw(4)

    this%kms = 1 - this%rwv      ; this%kds = 1
    this%kme = this%nz + this%rwv; this%kde = this%nz

    this%ipds = 1; this%ipde = this%nc
    this%jpds = 1; this%jpde = this%nc

    this%irs = max(this%ims, this%ids - this%rw - 1)
    this%ire = min(this%ime, this%ide + this%rw + 1)
    this%jrs = max(this%jms, this%jds - this%rw - 1)
    this%jre = min(this%jme, this%jde + this%rw + 1)
    
    this%its = merge(this%ids - 1, this%ids, .not. this%touch_edges(left)  )
    this%ite = merge(this%ide + 1, this%ide, .not. this%touch_edges(right) )
    this%jte = merge(this%jde + 1, this%jde, .not. this%touch_edges(top)   )
    this%jts = merge(this%jds - 1, this%jds, .not. this%touch_edges(bottom))

    this%igs = max(this%ipds, this%ids - 2 * this%rw - 1, this%ims)
    this%ige = min(this%ipde, this%ide + 2 * this%rw + 1, this%ime)
    this%jgs = max(this%jpds, this%jds - 2 * this%rw - 1, this%jms)
    this%jge = min(this%jpde, this%jde + 2 * this%rw + 1, this%jme)

    if (this%ngb(1)%ptr%active) this%ige = this%ide
    if (this%ngb(2)%ptr%active) this%igs = this%ids
    if (this%ngb(3)%ptr%active) this%jge = this%jde
    if (this%ngb(4)%ptr%active) this%jgs = this%jds

    allocate(this%r    (      this%npt               ,this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%x    (      this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%y    (      this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%z    (      this%npt               ,this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%eta  (      this%npt               ,                                    this%kms:this%kme)); this%eta = 0
    allocate(this%xi   (      this%npt               ,                                    this%kms:this%kme))
    allocate(this%lon  (      this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%lat  (      this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%Ah   (2,2,  this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%iAh  (2,2,  this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%Gh   (2,2,  this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%iGh  (2,2,  this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%Jh   (      this%npt               ,this%ims:this%ime,this%jms:this%jme                  ))
    allocate(this%CS   (3,3,3,this%pqs:this%pqe      ,this%ids:this%ide,this%jds:this%jde                  ))
    allocate(this%iGv  (2,3:3,this%npt               ,this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%Jv   (      this%npt               ,this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%Jab  (      this%npt               ,this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%iG   (3,3:3,this%pes(5):this%pee(6),this%ims:this%ime,this%jms:this%jme,this%kms:this%kme))
    allocate(this%omg  (3,    this%npt               ,this%ids:this%ide,this%jds:this%jde,this%kds:this%kde))

    allocate(this%ghost(                              this%ims:this%ime,this%jms:this%jme                  ))

    ! Set cell type.
    allocate(this%cell_type(this%ims:this%ime,this%jms:this%jme))

    this%cell_type = domain_cell
    if (this%touch_edges(left)) then
      this%cell_type(this%ims:this%ids-1,this%jms:this%jme) = left_halo_cell
    end if
    if (this%touch_edges(bottom)) then
      this%cell_type(this%ims:this%ime,this%jms:this%jds-1) = bottom_halo_cell
    end if
    if (this%touch_edges(right)) then
      this%cell_type(this%ide+1:this%ime,this%jms:this%jme) = right_halo_cell
    end if
    if (this%touch_edges(top)) then
      this%cell_type(this%ims:this%ime,this%jde+1:this%jme) = top_halo_cell
    end if
    if (this%touch_edges(left) .and. this%touch_edges(bottom)) then
      this%cell_type(this%ims:this%ids-1,this%jms:this%jds-1) = corner_cell
    end if
    if (this%touch_edges(bottom) .and. this%touch_edges(right)) then
      this%cell_type(this%ide+1:this%ime,this%jms:this%jds-1) = corner_cell
    end if
    if (this%touch_edges(right) .and. this%touch_edges(top)) then
      this%cell_type(this%ide+1:this%ime,this%jde+1:this%jme) = corner_cell
    end if
    if (this%touch_edges(top) .and. this%touch_edges(left)) then
      this%cell_type(this%ims:this%ids-1,this%jde+1:this%jme) = corner_cell
    end if
    
    select case (this%proj_type)
    case (equiangular)
      call this%init_equiangular()
    case default
      call log_notice('Unknown proj_type.', pid=this%pid)
      stop
    end select

    this%initialized = .true.

    ! Initialize neighbour domains
    if (this%domain_type == global_domain .and. this%active) then
      if (all(this%touch_edges)) then
        select case (this%ngb(1)%code)
        case (14, 21, 32, 43)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [hw,hw,rw,rw]
        case (64)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [rw,rw,hw,hw]
        case (54)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,rw,hw,hw]
        end select
        call this%ngb(1)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(1)%code)
        select case (this%ngb(2)%code)
        case (12, 23, 34, 41)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,hw,rw,rw]
        case (52)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,rw,hw,hw]
        case (62)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [rw,rw,hw,hw]
        end select
        call this%ngb(2)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(2)%code)
        select case (this%ngb(3)%code)
        case (16, 51)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,rw,hw,hw]
        case (26)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw,hw,rw,rw]
        case (46)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,hw,rw,rw]
        case (36, 63)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [rw,rw,hw,hw]
        end select
        call this%ngb(3)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(3)%code)
        select case (this%ngb(4)%code)
        case (15, 61)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [rw,rw,hw,hw]
        case (45)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,hw,rw,rw]
        case (25)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw,hw,rw,rw]
        case (35, 53)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,rw,hw,hw]
        end select
        call this%ngb(4)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(4)%code)
      ! Left-bottom corner
      else if (this%touch_edges(left) .and. this%touch_edges(bottom)) then
        ex_depth = min(max_ex_depth, merge(this%ngb(1)%ex_depth, max_ex_depth, this%ngb(1)%ex_depth > 0))
        select case (this%ngb(1)%code)
        case (14, 21, 32, 43)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [ex_depth,hw,rw,hw+ex_width]
        case (64)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [rw,hw+ex_width,hw,ex_depth]
        case (54)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,rw,ex_depth,hw]
        end select
        call this%ngb(1)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(1)%code)
        ex_depth = min(max_ex_depth, merge(this%ngb(3)%ex_depth, max_ex_depth, this%ngb(3)%ex_depth > 0))
        select case (this%ngb(3)%code)
        case (16, 51)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,hw+ex_width,ex_depth,hw]
        case (26)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [ex_depth,hw,hw+ex_width,rw]
        case (46)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,ex_depth,rw,hw+ex_width]
        case (36, 63)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,rw,hw,ex_depth]
        end select
        call this%ngb(3)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(3)%code)
      ! Bottom inner
      else if (.not. this%touch_edges(left) .and. .not. this%touch_edges(right) .and. this%touch_edges(bottom)) then
        select case (this%ngb(3)%code)
        case (16, 51)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = this%ids     ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        case (26)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-this%ide+1; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (46)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = this%ids     ; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (36, 63)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-this%ide+1; ngb_jds = 1            ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        end select
        call this%ngb(3)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(3)%code)
      ! Right-bottom corner
      else if (this%touch_edges(right) .and. this%touch_edges(bottom)) then
        ex_depth = min(max_ex_depth, merge(this%ngb(2)%ex_depth, max_ex_depth, this%ngb(2)%ex_depth > 0))
        select case (this%ngb(2)%code)
        case (12, 23, 34, 41)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,ex_depth,rw,hw+ex_width]
        case (52)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,hw+ex_width,ex_depth,hw]
        case (62)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,rw,hw,ex_depth]
        end select
        call this%ngb(2)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(2)%code)
        ex_depth = min(max_ex_depth, merge(this%ngb(3)%ex_depth, max_ex_depth, this%ngb(3)%ex_depth > 0))
        select case (this%ngb(3)%code)
        case (16, 51)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = this%ids     ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,rw,ex_depth,hw]
        case (26)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [ex_depth,hw,rw,hw+ex_width]
        case (46)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw,ex_depth,hw+ex_width,rw]
        case (36, 63)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [rw,hw+ex_width,hw,ex_depth]
        end select
        call this%ngb(3)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(3)%code)
      ! Right inner
      else if (this%touch_edges(right) .and. .not. this%touch_edges(bottom) .and. .not. this%touch_edges(top)) then
        select case (this%ngb(2)%code)
        case (12, 23, 34, 41)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = 1            ; ngb_jds = this%jds     ; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (52)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = this%jds     ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        case (62)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-this%jde+1; ngb_jds = 1            ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        end select
        call this%ngb(2)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(2)%code)
      ! Right-top corner
      else if (this%touch_edges(right) .and. this%touch_edges(top)) then
        ex_depth = min(max_ex_depth, merge(this%ngb(2)%ex_depth, max_ex_depth, this%ngb(2)%ex_depth > 0))
        select case (this%ngb(2)%code)
        case (12, 23, 34, 41)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw,ex_depth,hw+ex_width,rw]
        case (52)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,rw,ex_depth,hw]
        case (62)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [rw,hw+ex_width,hw,ex_depth]
        end select
        call this%ngb(2)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(2)%code)
        ex_depth = min(max_ex_depth, merge(this%ngb(4)%ex_depth, max_ex_depth, this%ngb(4)%ex_depth > 0))
        select case (this%ngb(4)%code)
        case (15, 61)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,rw,hw,ex_depth]
        case (45)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [hw,ex_depth,rw,hw+ex_width]
        case (25)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [ex_depth,hw,hw+ex_width,rw]
        case (35, 53)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,hw+ex_width,ex_depth,hw]
        end select
        call this%ngb(4)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(4)%code)
      ! Top inner
      else if (.not. this%touch_edges(left) .and. .not. this%touch_edges(right) .and. this%touch_edges(top)) then
        select case (this%ngb(4)%code)
        case (15, 61)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = this%ids     ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        case (45)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = nc-this%ide+1; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (25)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = this%ids     ; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (35, 53)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-this%ide+1; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        end select
        call this%ngb(4)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(4)%code)
      ! Left-top corner
      else if (this%touch_edges(left) .and. this%touch_edges(top)) then
        ex_depth = min(max_ex_depth, merge(this%ngb(1)%ex_depth, max_ex_depth, this%ngb(1)%ex_depth > 0))
        select case (this%ngb(1)%code)
        case (14, 21, 32, 43)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [ex_depth,hw,hw+ex_width,rw]
        case (64)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,rw,hw,ex_depth]
        case (54)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [rw,hw+ex_width,ex_depth,hw]
        end select
        call this%ngb(1)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(1)%code)
        ex_depth = min(max_ex_depth, merge(this%ngb(4)%ex_depth, max_ex_depth, this%ngb(4)%ex_depth > 0))
        select case (this%ngb(4)%code)
        case (15, 61)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = 1            ; ngb_jds = 1            ; ngb_hw = [rw,hw+ex_width,hw,ex_depth]
        case (45)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = 1            ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw,ex_depth,hw+ex_width,rw]
        case (25)
          ngb_nx = depth   ; ngb_ny = this%nbx; ngb_ids = nc-ngb_nx+1  ; ngb_jds = 1            ; ngb_hw = [ex_depth,hw,rw,hw+ex_width]
        case (35, 53)
          ngb_nx = this%nbx; ngb_ny = depth   ; ngb_ids = nc-ngb_nx+1  ; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,rw,ex_width,hw]
        end select
        call this%ngb(4)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(4)%code)
      ! Left inner
      else if (this%touch_edges(left) .and. .not. this%touch_edges(top) .and. .not. this%touch_edges(bottom)) then
        select case (this%ngb(1)%code)
        case (14, 21, 32, 43)
          ngb_nx = depth   ; ngb_ny = this%nby; ngb_ids = nc-ngb_nx+1  ; ngb_jds = this%jds     ; ngb_hw = [hw,hw,hw+ex_width,hw+ex_width]
        case (64)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = this%jds     ; ngb_jds = 1            ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        case (54)
          ngb_nx = this%nby; ngb_ny = depth   ; ngb_ids = nc-this%jde+1; ngb_jds = nc-ngb_ny+1  ; ngb_hw = [hw+ex_width,hw+ex_width,hw,hw]
        end select
        call this%ngb(1)%ptr%init(proj_type, hdim, vdim, nc, ndx, ndy, hws=ngb_hw, rw=this%rw, rwv=this%rwv, neq=neq, xeq=xeq, &
                                  ids=ngb_ids, jds=ngb_jds, nbx=ngb_nx, nby=ngb_ny, nz=this%nz, pid=pid, &
                                  ngb_code=this%ngb(1)%code)
      end if
    end if

  end subroutine cubed_sphere_domain_init

  subroutine cubed_sphere_domain_init_equiangular(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    integer i, j, k, p, q, m, n, idim
    integer count

    this%xmin = -pi0p25
    this%xmax =  pi0p25
    this%ymin = -pi0p25
    this%ymax =  pi0p25

    this%dx  = pi0p5 / this%nc
    this%dy  = pi0p5 / this%nc
    this%x0  = this%xmin + this%dx * (this%ids - 0.5_r8)
    this%y0  = this%ymin + this%dy * (this%jds - 0.5_r8)
    this%xds = this%x0 - this%dx * 0.5_r8
    this%yds = this%y0 - this%dy * 0.5_r8
    this%xde = this%xds + this%dx * this%nbx
    this%yde = this%yds + this%dy * this%nby
    
    this%deta = 1
    
    do j = this%jms, this%jme
      do i = this%ims, this%ime
        ! Cell center
        !   o---------------o
        !   |               |
        !   |               |
        !   |               |
        !   |       1       |
        !   |               |
        !   |               |
        !   |               |
        !   o---------------o
        this%x(this%pc,i,j) = this%x0 + (i - this%ids) * this%dx
        this%y(this%pc,i,j) = this%y0 + (j - this%jds) * this%dy
        
        ! Vertices
        ! 5 o---------------o 4
        !   |               |
        !   |               |
        !   |               |
        !   |               |
        !   |               |
        !   |               |
        !   |               |
        ! 2 o---------------o 3
        do idim = 1, this%hdim-1
          this%x((idim-1)*4+2,i,j) = this%x(this%pc,i,j) - 0.5d0 * this%dx
          this%y((idim-1)*4+2,i,j) = this%y(this%pc,i,j) - 0.5d0 * this%dy
          this%x((idim-1)*4+3,i,j) = this%x(this%pc,i,j) + 0.5d0 * this%dx
          this%y((idim-1)*4+3,i,j) = this%y(this%pc,i,j) - 0.5d0 * this%dy
          this%x((idim-1)*4+4,i,j) = this%x(this%pc,i,j) + 0.5d0 * this%dx
          this%y((idim-1)*4+4,i,j) = this%y(this%pc,i,j) + 0.5d0 * this%dy
          this%x((idim-1)*4+5,i,j) = this%x(this%pc,i,j) - 0.5d0 * this%dx
          this%y((idim-1)*4+5,i,j) = this%y(this%pc,i,j) + 0.5d0 * this%dy
        end do
        
        ! Edge quadrature with 2 points as an example
        !     10        11
        !   o--x---------x--o
        !   |               |
        !13 x               x 9
        !   |               |
        !   |               |
        !   |               |
        !12 x               x 8
        !   |               |
        !   o--x---------x--o
        !      6         7
        
        ! Left Face Points
        this%x(this%pes(1):this%pee(1),i,j) = this%x(this%pvs+0,i,j)
        do m = 1, this%neq**(this%hdim-2) 
          this%y(this%pes(1) + (m-1)*this%neq : this%pes(1) + m*this%neq-1, i, j) = this%y(this%pvs+0,i,j) + this%xeq * this%dy
        end do
        
        ! Right Face Points
        this%x(this%pes(2):this%pee(2),i,j) = this%x(this%pvs+1,i,j)
        do m = 1, this%neq**(this%hdim-2) 
          this%y(this%pes(2) + (m-1)*this%neq : this%pes(2) + m*this%neq-1, i,j) = this%y(this%pvs+1,i,j) + this%xeq * this%dy
        end do
        
        ! Bottom Face Points
        do m = 1, this%neq**(this%hdim-2) 
          this%x(this%pes(3) + (m-1)*this%neq : this%pes(3) + m*this%neq-1, i, j) = this%x(this%pvs+0,i,j) + this%xeq * this%dx
        end do
        this%y(this%pes(3):this%pee(3),i,j) = this%y(this%pvs+0,i,j)
        
        ! Top Face Points
        do m = 1, this%neq**(this%hdim-2) 
          this%x(this%pes(4) + (m-1)*this%neq : this%pes(4) + m*this%neq-1, i, j) = this%x(this%pvs+3,i,j) + this%xeq * this%dx
        end do
        this%y(this%pes(4):this%pee(4),i,j) = this%y(this%pvs+3,i,j)
        
        if(this%nfqv==1)then
          ! Down Face Points
          this%x(this%pes(5),i,j) = this%x(this%pc,i,j)
          this%y(this%pes(5),i,j) = this%y(this%pc,i,j)
          ! Up Face Points
          this%x(this%pes(6),i,j) = this%x(this%pc,i,j)
          this%y(this%pes(6),i,j) = this%y(this%pc,i,j)
        else
          ! Down Face Points
          k = this%pes(5)
          do q = 1, this%neq
            do p = 1, this%neq
              this%x(k,i,j) = this%x(2,i,j) + this%xeq(p) * this%dx
              this%y(k,i,j) = this%y(2,i,j) + this%xeq(q) * this%dy
              k = k + 1
            end do
          end do
          ! Up Face Points
          k = this%pes(6)
          do q = 1, this%neq
            do p = 1, this%neq
              this%x(k,i,j) = this%x(2,i,j) + this%xeq(p) * this%dx
              this%y(k,i,j) = this%y(2,i,j) + this%xeq(q) * this%dy
              k = k + 1
            end do
          end do
        endif
        
        ! Cell quadrature with 4 points as an example
        !   o---------------o
        !   | 16        17  |
        !   |  x         x  |
        !   |               |
        !   |               |
        !   |               |
        !   |  x         x  |
        !   | 14        15  |
        !   o---------------o
        k = this%pqs
        do m = 1, this%neq**(this%hdim-2)
          do q = 1, this%neq
            do p = 1, this%neq
              this%x(k,i,j) = this%x(2,i,j) + this%xeq(p) * this%dx
              this%y(k,i,j) = this%y(2,i,j) + this%xeq(q) * this%dy
              k = k + 1
            end do
          end do
        end do
        
        do p = 1, this%npt
          call this%cube_to_sphere(this%id, this%x(p,i,j), this%y(p,i,j), this%lon(p,i,j), this%lat(p,i,j))
        end do
      end do
    end do
    
    if(this%hdim==2)then
      do k = this%kms, this%kme
        ! Center
        this%eta(this%pc,k) = k - 0.5
        ! Vertices
        this%eta(this%pvs:this%pve,k) = this%eta(this%pc,k)
        ! Lateral
        do p = 1, 4
          this%eta(this%pes(p):this%pee(p),k) = this%eta(this%pc,k)
        end do
        ! Down
        this%eta(this%pes(5):this%pee(5),k) = k - 1
        ! Up
        this%eta(this%pes(6):this%pee(6),k) = k
        ! Quadrature
        this%eta(this%pqs:this%pqe,k) = this%eta(this%pc,k)
      enddo
    elseif(this%hdim==3)then
      do k = this%kms, this%kme
        ! Center
        this%eta(this%pc,k) = k - 0.5
        ! Vertices
        this%eta(this%pvs:this%pvs+3,k) = k - 1
        this%eta(this%pvs+4:this%pve,k) = k
        ! Lateral
        do m = 1, 4
          n = 0
          do q = 1, this%neq
            do p = 1, this%neq
              n = n + 1
              this%eta(this%pes(m)+n-1,k) = k - 1 + this%xeq(q)
            end do
          end do
        end do
        ! Down
        this%eta(this%pes(5):this%pee(5),k) = k - 1
        ! Up
        this%eta(this%pes(6):this%pee(6),k) = k
        ! Quadrature
        n = 0
        do q = 1, this%neq
          do p = 1, this%nfq
            n = n + 1
            this%eta(this%pqs+n-1,k) = k - 1 + this%xeq(q)
          end do
        end do
        ! Ghost points has to be set later
      enddo
    endif
    
    ! Set halo cells indices
    count = 0
    do j = this%jms, this%jme
      do i = this%ims, this%ime
        if( this%cell_type(i,j) /= domain_cell )then
          count = count + 1
        endif
      enddo
    enddo
    this%n_halo_cells = count
    
    allocate( this%halo_ij( 2, this%n_halo_cells ) )
    
    count = 0
    do j = this%jms, this%jme
      do i = this%ims, this%ime
        if( this%cell_type(i,j) /= domain_cell )then
          count = count + 1
          this%halo_ij(1,count) = i
          this%halo_ij(2,count) = j
        endif
      enddo
    enddo
    
  end subroutine cubed_sphere_domain_init_equiangular

  subroutine cubed_sphere_domain_set_horizontal_metrics(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    real(r16) :: r_r16
    real(r16) :: x_r16
    real(r16) :: y_r16
    real(r16) :: lon_r16
    real(r16) :: lat_r16
    real(r16) :: Gh_r16 (3,3)
    real(r16) :: iGh_r16(3,3)
    real(r16) :: Ah_r16 (3,3)
    real(r16) :: iAh_r16(3,3)
    real(r16) :: Jh_r16
    real(r16) :: CS_r16 (3,3,3)
    integer i, j, k, p
    
    call log_notice('Calculate horizontal metrics on cubed sphere.', pid=this%pid)

    !$OMP PARALLEL
    !$OMP DO PRIVATE(p,i,j,r_r16,x_r16,y_r16,lon_r16,lat_r16,Gh_r16,iGh_r16,Ah_r16,iAh_r16,Jh_r16,CS_r16) COLLAPSE(3)
    do j = this%jms, this%jme
      do i = this%ims, this%ime
        do p = 1, this%npt
          r_r16   = 1 !this%r  (p,i,j,k)
          x_r16   = this%x  (p,i,j)
          y_r16   = this%y  (p,i,j)
          lon_r16 = this%lon(p,i,j)
          lat_r16 = this%lat(p,i,j)
          call this%horizontal_metrics( this%id, r_r16,  &                
                                        x_r16,   y_r16,  &            
                                        lon_r16, lat_r16,&          
                                        Gh_r16,  iGh_r16,&
                                        Ah_r16,  iAh_r16,&
                                        Jh_r16,  CS_r16  )
          this%Gh (  :,:,p,i,j) = Gh_r16 (1:2,1:2)
          this%iGh(  :,:,p,i,j) = iGh_r16(1:2,1:2)
          this%Ah (  :,:,p,i,j) = Ah_r16 (1:2,1:2)
          this%iAh(  :,:,p,i,j) = iAh_r16(1:2,1:2)
          this%Jh (      p,i,j) = Jh_r16 
          if(   p>=this%pqs .and. p<=this%pqe &
          .and. i>=this%ids .and. i<=this%ide &
          .and. j>=this%jds .and. j<=this%jde &
          .and. k>=this%kds .and. k<=this%kde )then
            this%CS (:,:,:,p,i,j) = CS_r16
          endif
        end do
      end do
    end do
    !$OMP END DO NOWAIT

    !$OMP DO PRIVATE(p,i,j,r_r16,x_r16,y_r16,lon_r16,lat_r16,Gh_r16,iGh_r16,Ah_r16,iAh_r16,Jh_r16,CS_r16) COLLAPSE(2)
    do j = this%jgs, this%jge
      do i = this%igs, this%ige
        do p = 1, this%ghost(i,j)%ngp
          r_r16   = 1 !this%ghost(i,j)%r  (p,k)
          x_r16   = this%ghost(i,j)%x  (p)
          y_r16   = this%ghost(i,j)%y  (p)
          lon_r16 = this%ghost(i,j)%lon(p)
          lat_r16 = this%ghost(i,j)%lat(p)
          call this%horizontal_metrics( this%id, r_r16,  &                
                                        x_r16,   y_r16,  &            
                                        lon_r16, lat_r16,&          
                                        Gh_r16,  iGh_r16,&
                                        Ah_r16,  iAh_r16,&
                                        Jh_r16,  CS_r16  )
          this%ghost(i,j)%Gh (:,:,p) = Gh_r16 (1:2,1:2)
          this%ghost(i,j)%iGh(:,:,p) = iGh_r16(1:2,1:2)
          this%ghost(i,j)%Ah (:,:,p) = Ah_r16 (1:2,1:2)
          this%ghost(i,j)%iAh(:,:,p) = iAh_r16(1:2,1:2)
          this%ghost(i,j)%Jh (    p) = Jh_r16 
        end do
      end do
    end do
    !$OMP END DO
    !$OMP END PARALLEL
  end subroutine cubed_sphere_domain_set_horizontal_metrics

  subroutine cubed_sphere_domain_set_vertical_metrics(this,Ra,Rb,Rr,p,i,j,k)

    class(cubed_sphere_domain_type), intent(inout) :: this
    real(r8) :: Ra, Rb, Rr
    integer p, i, j, k
    
    this%iGv(1,3,p,i,j,k) = -Ra / Rr
    this%iGv(2,3,p,i,j,k) = -Rb / Rr
    this%Jv (    p,i,j,k) = Rr
    
  end subroutine cubed_sphere_domain_set_vertical_metrics

  subroutine cubed_sphere_domain_set_ghost_vertical_metrics(this,Ra,Rb,Rr,p,i,j,k)

    class(cubed_sphere_domain_type), intent(inout) :: this
    real(r8) :: Ra, Rb, Rr
    integer p, i, j, k
    
    this%ghost(i,j)%Jv(p,k) = Rr
    
  end subroutine cubed_sphere_domain_set_ghost_vertical_metrics

  subroutine cubed_sphere_domain_count_ghost_points(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    integer is, ie, js, je, ks, ke
    integer i, j, k, p, it, jt
    integer idom1, idom2, idom
    real(r8) x1, y1, x2, y2, wgt

    !if (.not. this%active) then
    !  if (this%ngb(left)%ptr%active) then
    !    is = this%irs; ie = this%ide
    !    js = this%jrs; je = this%jre
    !  else if (this%ngb(right)%ptr%active) then
    !    is = this%ids; ie = this%ire
    !    js = this%jrs; je = this%jre
    !  else if (this%ngb(bottom)%ptr%active) then
    !    is = this%irs; ie = this%ire
    !    js = this%jrs; je = this%jde
    !  else if (this%ngb(top)%ptr%active) then
    !    is = this%irs; ie = this%ire
    !    js = this%jds; je = this%jre
    !  end if
    !else
    !  is = this%irs; ie = this%ire
    !  js = this%jrs; je = this%jre
    !end if
    !ks = this%kms; ke = this%kme

    if (.not. this%active) then
      if (this%ngb(left)%ptr%active) then
        is = this%ims; ie = this%ide
        js = this%jms; je = this%jme
      else if (this%ngb(right)%ptr%active) then
        is = this%ids; ie = this%ime
        js = this%jms; je = this%jme
      else if (this%ngb(bottom)%ptr%active) then
        is = this%ims; ie = this%ime
        js = this%jms; je = this%jde
      else if (this%ngb(top)%ptr%active) then
        is = this%ims; ie = this%ime
        js = this%jds; je = this%jme
      end if
    else
      is = this%ims; ie = this%ime
      js = this%jms; je = this%jme
    end if
    ks = this%kms; ke = this%kme
    
    ! Count ghost point number
    do j = js, je
      do i = is, ie
        select case (this%cell_type(i,j))
        case (left_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(1)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
          end do
        case (right_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(2)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
          enddo
        case (bottom_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(3)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
          end do
        case (top_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(4)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
          end do
        case (corner_cell)
          do p = this%pqs, this%pqe
            call this%locate_domain(this%lon(p,i,j), this%lat(p,i,j), idom1, x1, y1, idom2, x2, y2)
            if (merge(this%ngb(1)%ptr%id == idom1, .false., associated(this%ngb(1)%ptr))) then
              call this%ngb(1)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(3)%ptr%id, .false., associated(this%ngb(3)%ptr))) then
                call this%ngb(3)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(4)%ptr%id, .false., associated(this%ngb(4)%ptr))) then
                call this%ngb(4)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              end if
            else if (merge(this%ngb(2)%ptr%id == idom1, .false., associated(this%ngb(2)%ptr))) then
              call this%ngb(2)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(3)%ptr%id, .false., associated(this%ngb(3)%ptr))) then
                call this%ngb(3)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(4)%ptr%id, .false., associated(this%ngb(4)%ptr))) then
                call this%ngb(4)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              end if
            else if (merge(this%ngb(3)%ptr%id == idom1, .false., associated(this%ngb(3)%ptr))) then
              call this%ngb(3)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(1)%ptr%id, .false., associated(this%ngb(1)%ptr))) then
                call this%ngb(1)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(2)%ptr%id, .false., associated(this%ngb(2)%ptr))) then
                call this%ngb(2)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              end if
            else if (merge(this%ngb(4)%ptr%id == idom1, .false., associated(this%ngb(4)%ptr))) then
              call this%ngb(4)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(1)%ptr%id, .false., associated(this%ngb(1)%ptr))) then
                call this%ngb(1)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(2)%ptr%id, .false., associated(this%ngb(2)%ptr))) then
                call this%ngb(2)%ptr%count_ghost_point(this%lon(p,i,j), this%lat(p,i,j))
              end if
            end if
          end do
        end select
      end do
    end do
    
  end subroutine cubed_sphere_domain_count_ghost_points

  subroutine cubed_sphere_domain_allocate_ghost_points(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    integer it, jt
    
    do jt = this%jms, this%jme
      do it = this%ims, this%ime
        allocate( this%ghost(it,jt)%idom        (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%ps          (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%is          (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%js          (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%wgt         (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%x           (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%y           (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%lon         (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%lat         (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%r           (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        allocate( this%ghost(it,jt)%z           (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        allocate( this%ghost(it,jt)%eta         (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        allocate( this%ghost(it,jt)%xi          (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        allocate( this%ghost(it,jt)%Ah          (2,2,  this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%iAh         (2,2,  this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%Gh          (2,2,  this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%iGh         (2,2,  this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%Jh          (      this%ghost(it,jt)%ngp                  ) )
        allocate( this%ghost(it,jt)%Jv          (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        allocate( this%ghost(it,jt)%Jab         (      this%ghost(it,jt)%ngp,this%kms:this%kme) )
        
        ! Reset ngp for add ghost point info
        this%ghost(it,jt)%ngp = 0
      enddo
    enddo

  end subroutine cubed_sphere_domain_allocate_ghost_points

  subroutine cubed_sphere_domain_add_ghost_points(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    integer is, ie, js, je, ks, ke
    integer i, j, k, p, it, jt
    integer idom1, idom2, idom
    real(r8) x1, y1, x2, y2, wgt

    !if (.not. this%active) then
    !  if (this%ngb(left)%ptr%active) then
    !    is = this%irs; ie = this%ide
    !    js = this%jrs; je = this%jre
    !  else if (this%ngb(right)%ptr%active) then
    !    is = this%ids; ie = this%ire
    !    js = this%jrs; je = this%jre
    !  else if (this%ngb(bottom)%ptr%active) then
    !    is = this%irs; ie = this%ire
    !    js = this%jrs; je = this%jde
    !  else if (this%ngb(top)%ptr%active) then
    !    is = this%irs; ie = this%ire
    !    js = this%jds; je = this%jre
    !  end if
    !else
    !  is = this%irs; ie = this%ire
    !  js = this%jrs; je = this%jre
    !end if
    !ks = this%kms; ke = this%kme

    if (.not. this%active) then
      if (this%ngb(left)%ptr%active) then
        is = this%ims; ie = this%ide
        js = this%jms; je = this%jme
      else if (this%ngb(right)%ptr%active) then
        is = this%ids; ie = this%ime
        js = this%jms; je = this%jme
      else if (this%ngb(bottom)%ptr%active) then
        is = this%ims; ie = this%ime
        js = this%jms; je = this%jde
      else if (this%ngb(top)%ptr%active) then
        is = this%ims; ie = this%ime
        js = this%jds; je = this%jme
      end if
    else
      is = this%ims; ie = this%ime
      js = this%jms; je = this%jme
    end if
    ks = this%kms; ke = this%kme
    
    ! Add point info
    do j = js, je
      do i = is, ie
        select case (this%cell_type(i,j))
        case (left_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(1)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), 1._r8, it, jt)
            if (.not. this%active) then
              select case (10 * this%ngb_code + this%ngb(1)%ptr%id)
              case (154, 164)
                this%ngb(1)%ex_depth = max(this%ngb(1)%ex_depth, this%ngb(1)%ptr%ide - it + 1)
              case (354, 364)
                this%ngb(1)%ex_depth = max(this%ngb(1)%ex_depth,                       it    )
              case (543, 514, 521, 532)
                this%ngb(1)%ex_depth = max(this%ngb(1)%ex_depth, this%ngb(1)%ptr%jde - jt + 1)
              case (614, 643, 632, 621)
                this%ngb(1)%ex_depth = max(this%ngb(1)%ex_depth,                       jt    )
              end select
            end if
          end do
        case (right_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(2)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), 1._r8, it, jt)
            if (.not. this%active) then
              select case (10 * this%ngb_code + this%ngb(2)%ptr%id)
              case (152, 162)
                this%ngb(2)%ex_depth = max(this%ngb(2)%ex_depth,                       it    )
              case (352, 362)
                this%ngb(2)%ex_depth = max(this%ngb(2)%ex_depth, this%ngb(2)%ptr%ide - it + 1)
              case (534, 541, 512, 523)
                this%ngb(2)%ex_depth = max(this%ngb(2)%ex_depth, this%ngb(2)%ptr%jde - jt + 1)
              case (612, 641, 634, 623)
                this%ngb(2)%ex_depth = max(this%ngb(2)%ex_depth,                       jt    )
              end select
            end if
          end do
        case (bottom_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(3)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), 1._r8, it, jt)
            if (.not. this%active) then
              select case (10 * this%ngb_code + this%ngb(3)%ptr%id)
              case (146, 126)
                this%ngb(3)%ex_depth = max(this%ngb(3)%ex_depth, this%ngb(3)%ptr%jde - jt + 1)
              case (251, 216, 236, 463)
                this%ngb(3)%ex_depth = max(this%ngb(3)%ex_depth, this%ngb(3)%ptr%ide - it + 1)
              case (263, 451, 436, 416)
                this%ngb(3)%ex_depth = max(this%ngb(3)%ex_depth,                       it    )
              case (326, 346)
                this%ngb(3)%ex_depth = max(this%ngb(3)%ex_depth,                       jt    )
              end select
            end if
          end do
        case (top_halo_cell)
          do p = this%pqs, this%pqe
            call this%ngb(4)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), 1._r8, it, jt)
            if (.not. this%active) then
              select case (10 * this%ngb_code + this%ngb(4)%ptr%id)
              case (125, 145)
                this%ngb(4)%ex_depth = max(this%ngb(4)%ex_depth,                       jt    )
              case (253, 435, 461, 415)
                this%ngb(4)%ex_depth = max(this%ngb(4)%ex_depth,                       it    )
              case (215, 261, 235, 453)
                this%ngb(4)%ex_depth = max(this%ngb(4)%ex_depth, this%ngb(4)%ptr%ide - it + 1)
              case (325, 345)
                this%ngb(4)%ex_depth = max(this%ngb(4)%ex_depth, this%ngb(4)%ptr%jde - jt + 1)
              end select
            end if
          end do
        case (corner_cell)
          do p = this%pqs, this%pqe
            call this%locate_domain(this%lon(p,i,j), this%lat(p,i,j), idom1, x1, y1, idom2, x2, y2)
            wgt = merge(0.5_r8, 1._r8, idom1 /= idom2)
            if (merge(this%ngb(1)%ptr%id == idom1, .false., associated(this%ngb(1)%ptr))) then
              call this%ngb(1)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(3)%ptr%id, .false., associated(this%ngb(3)%ptr))) then
                call this%ngb(3)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(4)%ptr%id, .false., associated(this%ngb(4)%ptr))) then
                call this%ngb(4)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              end if
            else if (merge(this%ngb(2)%ptr%id == idom1, .false., associated(this%ngb(2)%ptr))) then
              call this%ngb(2)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(3)%ptr%id, .false., associated(this%ngb(3)%ptr))) then
                call this%ngb(3)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(4)%ptr%id, .false., associated(this%ngb(4)%ptr))) then
                call this%ngb(4)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              end if
            else if (merge(this%ngb(3)%ptr%id == idom1, .false., associated(this%ngb(3)%ptr))) then
              call this%ngb(3)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(1)%ptr%id, .false., associated(this%ngb(1)%ptr))) then
                call this%ngb(1)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(2)%ptr%id, .false., associated(this%ngb(2)%ptr))) then
                call this%ngb(2)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              end if
            else if (merge(this%ngb(4)%ptr%id == idom1, .false., associated(this%ngb(4)%ptr))) then
              call this%ngb(4)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              if (idom1 /= idom2 .and. merge(idom2 == this%ngb(1)%ptr%id, .false., associated(this%ngb(1)%ptr))) then
                call this%ngb(1)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              else if (idom1 /= idom2 .and. merge(idom2 == this%ngb(2)%ptr%id, .false., associated(this%ngb(2)%ptr))) then
                call this%ngb(2)%ptr%add_ghost_point(this%id, p, i, j, k, this%lon(p,i,j), this%lat(p,i,j), this%eta(p,:), wgt, it, jt)
              end if
            end if
          end do
        end select
      end do
    end do
    
    if (.not. this%active) then
      if (this%ngb(1)%ex_depth > 0) this%ngb(1)%ex_depth = this%ngb(1)%ex_depth - this%rw - 1
      if (this%ngb(2)%ex_depth > 0) this%ngb(2)%ex_depth = this%ngb(2)%ex_depth - this%rw - 1
      if (this%ngb(3)%ex_depth > 0) this%ngb(3)%ex_depth = this%ngb(3)%ex_depth - this%rw - 1
      if (this%ngb(4)%ex_depth > 0) this%ngb(4)%ex_depth = this%ngb(4)%ex_depth - this%rw - 1
    end if

  end subroutine cubed_sphere_domain_add_ghost_points

  subroutine cubed_sphere_domain_set_ghost_ij(this)
    class(cubed_sphere_domain_type), intent(inout) :: this
    
    integer count
    integer i,j
    
    count = 0
    do j = this%jgs, this%jge
      do i = this%igs, this%ige
        if( this%ghost(i,j)%ngp > 0 )then
          count = count + 1
        end if
      end do
    end do
    this%n_gst_cells = count
    
    allocate(this%gst_ij(2,this%n_gst_cells))

    count = 0
    do j = this%jgs, this%jge
      do i = this%igs, this%ige
        if( this%ghost(i,j)%ngp > 0 )then
          count = count + 1
          this%gst_ij(:,count) = [i,j]
        end if
      end do
    end do
    
  end subroutine cubed_sphere_domain_set_ghost_ij
  
  subroutine cubed_sphere_domain_count_ghost_point(this, lon, lat)

    class(cubed_sphere_domain_type), intent(inout) :: this
    real(r8), intent(in) :: lon   ! Point longitude
    real(r8), intent(in) :: lat   ! Point latitude
    
    integer :: it   ! Target cell x-axis index
    integer :: jt   ! Target cell y-axis index
    
    real(r8) x, y

    call this%sphere_to_cube(this%id, lon, lat, x, y)
    it = min(max(floor((x - this%xmin) / this%dx) + 1, 1), this%ndx)
    jt = min(max(floor((y - this%ymin) / this%dy) + 1, 1), this%ndy)
    if (it < this%ims .or. it > this%ime .or. jt < this%jms .or. jt > this%jme) then
      call log_error('Ghost point is outside neighbor domain!', __FILE__, __LINE__)
    end if
    this%ghost(it,jt)%ngp = this%ghost(it,jt)%ngp + 1
    
  end subroutine cubed_sphere_domain_count_ghost_point

  subroutine cubed_sphere_domain_add_ghost_point(this, idom, ps, is, js, ks, lon, lat, eta, wgt, it, jt)

    class(cubed_sphere_domain_type), intent(inout) :: this
    integer , intent(in) :: idom  ! Source domain ID
    integer , intent(in) :: ps    ! Source point index
    integer , intent(in) :: is    ! Source cell x-axis index
    integer , intent(in) :: js    ! Source cell y-axis index
    integer , intent(in) :: ks    ! Source cell k-axis index
    real(r8), intent(in) :: lon   ! Point longitude
    real(r8), intent(in) :: lat   ! Point latitude
    real(r8), intent(in) :: eta(:)! Point eta
    real(r8), intent(in) :: wgt
    integer , intent(out) :: it   ! Target cell x-axis index
    integer , intent(out) :: jt   ! Target cell y-axis index
    
    real(r8) x, y
    integer p

    call this%sphere_to_cube(this%id, lon, lat, x, y)
    it = min(max(floor((x - this%xmin) / this%dx) + 1, 1), this%ndx)
    jt = min(max(floor((y - this%ymin) / this%dy) + 1, 1), this%ndy)
    if (it < this%ims .or. it > this%ime .or. jt < this%jms .or. jt > this%jme) then
      call log_error('Ghost point is outside neighbor domain!', __FILE__, __LINE__)
    end if
    this%ghost(it,jt)%ngp = this%ghost(it,jt)%ngp + 1
    if (this%ghost(it,jt)%ngp > this%max_ngp) then
      call log_error('Not sufficient array to store ghost points! Enlarge it!', __FILE__, __LINE__)
    end if
    p = this%ghost(it,jt)%ngp
    this%ghost(it,jt)%idom(p) = idom
    this%ghost(it,jt)%ps  (p) = ps
    this%ghost(it,jt)%is  (p) = is
    this%ghost(it,jt)%js  (p) = js
    this%ghost(it,jt)%wgt (p) = wgt
    this%ghost(it,jt)%lon (p) = lon
    this%ghost(it,jt)%lat (p) = lat
    this%ghost(it,jt)%x   (p) = x
    this%ghost(it,jt)%y   (p) = y
    this%ghost(it,jt)%eta (p,:) = eta

  end subroutine cubed_sphere_domain_add_ghost_point

  real(r8) function cubed_sphere_domain_edge_quad(this, x) result(res)

    class(cubed_sphere_domain_type), intent(in) :: this
    real(r8), intent(in) :: x(:)

    res = dot_product(this%weq_edge, x)

  end function cubed_sphere_domain_edge_quad

  real(r8) function cubed_sphere_domain_cell_quad(this, x) result(res)

    class(cubed_sphere_domain_type), intent(in) :: this
    real(r8), intent(in) :: x(:)

    res = dot_product(this%weq_cell, x)

  end function cubed_sphere_domain_cell_quad

  real(r8) function cubed_sphere_domain_edge_quad_ver(this, x) result(res)

    class(cubed_sphere_domain_type), intent(in) :: this
    real(r8), intent(in) :: x(:)

    res = dot_product(this%weq_v, x)

  end function cubed_sphere_domain_edge_quad_ver

  logical function cubed_sphere_domain_is_xy_in_domain(this, x, y) result(res)

    class(cubed_sphere_domain_type), intent(in) :: this
    real(r8), intent(in) :: x
    real(r8), intent(in) :: y

    res = .true.
    if (x < this%xds .or. x > this%xde .or. y < this%yds .or. y > this%yde) res = .false.

  end function cubed_sphere_domain_is_xy_in_domain

  subroutine cubed_sphere_domain_clear(this)

    class(cubed_sphere_domain_type), intent(inout) :: this

    if (allocated(this%eta     )) deallocate(this%eta     )
    if (allocated(this%xi      )) deallocate(this%xi      )
    if (allocated(this%r       )) deallocate(this%r       )
    if (allocated(this%xeq     )) deallocate(this%xeq     )
    if (allocated(this%weq1d   )) deallocate(this%weq1d   )
    if (allocated(this%weq2d   )) deallocate(this%weq2d   )
    if (allocated(this%weq3d   )) deallocate(this%weq3d   )
    if (allocated(this%weq_v   )) deallocate(this%weq_v   )
    if (allocated(this%weq_edge)) deallocate(this%weq_edge)
    if (allocated(this%weq_cell)) deallocate(this%weq_cell)
    if (allocated(this%x       )) deallocate(this%x       )
    if (allocated(this%y       )) deallocate(this%y       )
    if (allocated(this%z       )) deallocate(this%z       )
    if (allocated(this%lon     )) deallocate(this%lon     )
    if (allocated(this%lat     )) deallocate(this%lat     )
          
    if (allocated(this%Ah      )) deallocate(this%Ah      )
    if (allocated(this%iAh     )) deallocate(this%iAh     )
    if (allocated(this%Gh      )) deallocate(this%Gh      )
    if (allocated(this%iGh     )) deallocate(this%iGh     )
    if (allocated(this%Jh      )) deallocate(this%Jh      )
    if (allocated(this%CS      )) deallocate(this%CS      )
          
    if (allocated(this%iGv     )) deallocate(this%iGv     )
    if (allocated(this%Jv      )) deallocate(this%Jv      )
          
    if (allocated(this%Jab     )) deallocate(this%Jab     )
    if (allocated(this%iG      )) deallocate(this%iG      )
          
    if (allocated(this%omg     )) deallocate(this%omg     )
          
    if (allocated(this%ghost   )) deallocate(this%ghost   )

    if (allocated(this%cell_type)) deallocate(this%cell_type)

  end subroutine cubed_sphere_domain_clear

  subroutine cubed_sphere_domain_final(this)

    type(cubed_sphere_domain_type), intent(inout) :: this

    call this%clear()

  end subroutine cubed_sphere_domain_final

  subroutine ghost_clear(this)
    
    class(ghost_point_type), intent(inout) :: this

    if(allocated(this%idom        )) deallocate(this%idom        )
    if(allocated(this%ps          )) deallocate(this%ps          )
    if(allocated(this%is          )) deallocate(this%is          )
    if(allocated(this%js          )) deallocate(this%js          )
    if(allocated(this%wgt         )) deallocate(this%wgt         )
    if(allocated(this%x           )) deallocate(this%x           )
    if(allocated(this%y           )) deallocate(this%y           )
    if(allocated(this%lon         )) deallocate(this%lon         )
    if(allocated(this%lat         )) deallocate(this%lat         )
    if(allocated(this%r           )) deallocate(this%r           )
    if(allocated(this%z           )) deallocate(this%z           )
    if(allocated(this%eta         )) deallocate(this%eta         )
    if(allocated(this%xi          )) deallocate(this%xi          )
    if(allocated(this%Ah          )) deallocate(this%Ah          )
    if(allocated(this%iAh         )) deallocate(this%iAh         )
    if(allocated(this%Gh          )) deallocate(this%Gh          )
    if(allocated(this%iGh         )) deallocate(this%iGh         )
    if(allocated(this%Jh          )) deallocate(this%Jh          )
    if(allocated(this%Jv          )) deallocate(this%Jv          )
    if(allocated(this%Jab         )) deallocate(this%Jab         )
  end subroutine ghost_clear

  subroutine ghost_final(this)

    type(ghost_point_type), intent(inout) :: this
    
    call this%clear()

  end subroutine ghost_final

end module cubed_sphere_domain_mod
